On the importance of Dualization
using DynamicPolynomials
using SumOfSquaresSum-of-Squares programs are usually solved by SemiDefinite Programming solvers (SDPs). These programs can be represented into two different formats: Either the standard conic form, also known as kernel form:
\[\begin{aligned} \min\limits_{Q \in \mathbb{S}_n} & \langle C, Q \rangle\\ \text{subject to:} & \langle A_i, Q \rangle = b_i, \quad i=1,2,\ldots,m\\ & Q \succeq 0, \end{aligned}\]
or the geometric conic form, also known as image form:
\[\begin{aligned} \max\limits_{y \in \mathbb{R}^m} & \langle b, y \rangle\\ \text{subject to:} & C \succeq \sum_{i=1}^m A_i y_i\\ & y\ \mathsf{free}, \end{aligned}\]
In this tutorial, we investigate in which of these two forms a Sum-of-Squares constraint should be written into. Consider the simple example of trying to determine whether the following univariate polynomial is a Sum-of-Squares:
import SCS
@polyvar x
p = (x + 1)^2 * (x + 2)^2
model_scs = Model(SCS.Optimizer)
con_ref = @constraint(model_scs, p in SOSCone())
optimize!(model_scs)Success: SDP solved
Primal objective value: 0.0000000e+00
Dual objective value: 0.0000000e+00
Relative primal infeasibility: 4.15e-17
Relative dual infeasibility: 5.00e-11
Real Relative Gap: 0.00e+00
XZ Relative Gap: 3.63e-10
DIMACS error measures: 6.34e-17 0.00e+00 5.00e-11 0.00e+00 0.00e+00 3.63e-10
------------------------------------------------------------------
SCS v3.2.4 - Splitting Conic Solver
(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem: variables n: 6, constraints m: 11
cones: z: primal zero / dual free vars: 5
s: psd vars: 6, ssize: 1
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07
alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
max_iters: 100000, normalize: 1, rho_x: 1.00e-06
acceleration_lookback: 10, acceleration_interval: 10
lin-sys: sparse-direct-amd-qdldl
nnz(A): 12, nnz(P): 0
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 1.30e+01 1.87e+00 2.95e+01 1.47e+01 1.00e-01 7.85e-05
50| 6.01e-07 2.84e-07 1.95e-06 -9.76e-07 1.00e-01 1.82e-04
------------------------------------------------------------------
status: solved
timings: total: 1.83e-04s = setup: 4.75e-05s + solve: 1.35e-04s
lin-sys: 1.23e-05s, cones: 9.41e-05s, accel: 2.47e-06s
------------------------------------------------------------------
objective = -0.000001
------------------------------------------------------------------As we can see in the log, SCS reports 6 variables and 11 constraints. We can also choose to dualize the problem before it is passed to SCS as follows:
using Dualization
model_dual_scs = Model(dual_optimizer(SCS.Optimizer))
@objective(model_dual_scs, Max, 0.0)
con_ref = @constraint(model_dual_scs, p in SOSCone())
optimize!(model_dual_scs)------------------------------------------------------------------
SCS v3.2.4 - Splitting Conic Solver
(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem: variables n: 5, constraints m: 6
cones: s: psd vars: 6, ssize: 1
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07
alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
max_iters: 100000, normalize: 1, rho_x: 1.00e-06
acceleration_lookback: 10, acceleration_interval: 10
lin-sys: sparse-direct-amd-qdldl
nnz(A): 6, nnz(P): 0
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 1.97e+02 1.05e+01 4.39e+03 -2.20e+03 1.00e-01 7.11e-05
50| 6.29e-06 1.81e-06 2.22e-05 -1.11e-05 1.00e-01 1.69e-04
------------------------------------------------------------------
status: solved
timings: total: 1.70e-04s = setup: 4.08e-05s + solve: 1.29e-04s
lin-sys: 9.40e-06s, cones: 9.28e-05s, accel: 2.55e-06s
------------------------------------------------------------------
objective = -0.000011
------------------------------------------------------------------This time, SCS reports 5 variables and 6 constraints.
Bridges operating behind the scenes
The difference comes from the fact that, when designing the JuMP interface of SCS, it was decided that the model would be read in the image form. SCS therefore declares that it only supports free variables, represented in JuMP as variables in MOI.Reals and affine semidefinite constraints, represented in JuMP as MOI.VectorAffineFunction-in-MOI.PositiveSemidefiniteConeTriangle constraints. On the other hand, SumOfSquares gave the model in kernel form so the positive semidefinite (PSD) variables were reformulated as free variables constrained to be PSD using an affine PSD constraints.
This transformation is done transparently without warning but it can be inspected using print_active_bridges. As shown below, we can see Unsupported variable: MOI.PositiveSemidefiniteConeTriangle and adding as constraint indicating that PSD variables are not supported and they are added as free variables. Then we have Unsupported constraint: MOI.VectorOfVariables-in-MOI.PositiveSemidefiniteConeTriangle indicating that SCS does not support constraining variables in the PSD cone so it will just convert it into affine expressions in the PSD cone. Of course, this is equivalent but it means that SCS will not exploit this particular structure of the problem hence solving might be less efficient.
print_active_bridges(model_scs) * Supported objective: MOI.ScalarAffineFunction{Float64}
* Unsupported constraint: MOI.VectorAffineFunction{Float64}-in-SumOfSquares.SOSPolynomialSet{FullSpace, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, SumOfSquares.Certificate.NewtonFilter{SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}}}}}
| bridged by:
| SumOfSquares.Bridges.Constraint.SOSPolynomialBridge{Float64, MOI.VectorAffineFunction{Float64}, FullSpace, Union{MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.Nonnegatives}, MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle}, MOI.ConstraintIndex{MOI.VectorOfVariables, SumOfSquares.EmptyCone}, MOI.ConstraintIndex{MOI.VectorOfVariables, SumOfSquares.PositiveSemidefinite2x2ConeTriangle}, Vector{Union{MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.Nonnegatives}, MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle}, MOI.ConstraintIndex{MOI.VectorOfVariables, SumOfSquares.EmptyCone}, MOI.ConstraintIndex{MOI.VectorOfVariables, SumOfSquares.PositiveSemidefinite2x2ConeTriangle}}}}, Union{SumOfSquares.EmptyCone, SumOfSquares.PositiveSemidefinite2x2ConeTriangle, MOI.Nonnegatives, MOI.PositiveSemidefiniteConeTriangle}, MOI.PositiveSemidefiniteConeTriangle, MonomialBasis, MonomialBasis, SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, SumOfSquares.Certificate.NewtonFilter{SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}}}}, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}
| may introduce:
| * Unsupported constraint: MOI.VectorAffineFunction{Float64}-in-PolyJuMP.ZeroPolynomialSet{FullSpace, MonomialBasis, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}
| | bridged by:
| | PolyJuMP.Bridges.Constraint.ZeroPolynomialBridge{Float64, MOI.VectorAffineFunction{Float64}, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}
| | may introduce:
| | * Supported constraint: MOI.VectorAffineFunction{Float64}-in-MOI.Zeros
| * Unsupported variable: SumOfSquares.EmptyCone
| | adding as constraint:
| | * Supported variable: MOI.Reals
| | * Unsupported constraint: MOI.VectorOfVariables-in-SumOfSquares.EmptyCone
| | | bridged by:
| | | SumOfSquares.Bridges.Constraint.EmptyBridge{Float64, MOI.VectorOfVariables}
| | | may introduce:
| * Unsupported variable: MOI.Nonnegatives
| | adding as constraint:
| | * Supported variable: MOI.Reals
| | * Unsupported constraint: MOI.VectorOfVariables-in-MOI.Nonnegatives
| | | bridged by:
| | | MOIB.Constraint.FunctionConversionBridge{Float64, MOI.VectorAffineFunction{Float64}, MOI.VectorOfVariables, MOI.Nonnegatives}
| | | may introduce:
| | | * Supported constraint: MOI.VectorAffineFunction{Float64}-in-MOI.Nonnegatives
| * Unsupported variable: SumOfSquares.PositiveSemidefinite2x2ConeTriangle
| | bridged by:
| | SumOfSquares.Bridges.Variable.PositiveSemidefinite2x2Bridge{Float64}
| | may introduce:
| | * Unsupported variable: MOI.RotatedSecondOrderCone
| | | adding as constraint:
| | | * Supported variable: MOI.Reals
| | | * Unsupported constraint: MOI.VectorOfVariables-in-MOI.RotatedSecondOrderCone
| | | | bridged by:
| | | | MOIB.Constraint.RSOCtoSOCBridge{Float64, MOI.VectorAffineFunction{Float64}, MOI.VectorOfVariables}
| | | | may introduce:
| | | | * Supported constraint: MOI.VectorAffineFunction{Float64}-in-MOI.SecondOrderCone
| * Unsupported variable: MOI.PositiveSemidefiniteConeTriangle
| | adding as constraint:
| | * Supported variable: MOI.Reals
| | * Unsupported constraint: MOI.VectorOfVariables-in-MOI.PositiveSemidefiniteConeTriangle
| | | bridged by:
| | | MOIB.Constraint.SetDotScalingBridge{Float64, MOI.PositiveSemidefiniteConeTriangle, MOI.VectorAffineFunction{Float64}, MOI.VectorOfVariables}
| | | may introduce:
| | | * Unsupported constraint: MOI.VectorAffineFunction{Float64}-in-MOI.Scaled{MOI.PositiveSemidefiniteConeTriangle}
| | | | bridged by:
| | | | SCS.ScaledPSDConeBridge{Float64, MOI.VectorAffineFunction{Float64}}
| | | | may introduce:
| | | | * Supported constraint: MOI.VectorAffineFunction{Float64}-in-SCS.ScaledPSDConeWith the dual version, we can see that variables in the PSD cone are supported directly hence we don't need that extra conversion.
print_active_bridges(model_dual_scs) * Supported objective: MOI.ScalarAffineFunction{Float64}
* Unsupported constraint: MOI.VectorAffineFunction{Float64}-in-SumOfSquares.SOSPolynomialSet{FullSpace, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, SumOfSquares.Certificate.NewtonFilter{SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}}}}}
| bridged by:
| SumOfSquares.Bridges.Constraint.SOSPolynomialBridge{Float64, MOI.VectorAffineFunction{Float64}, FullSpace, Union{MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.Nonnegatives}, MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle}, MOI.ConstraintIndex{MOI.VectorOfVariables, SumOfSquares.EmptyCone}, MOI.ConstraintIndex{MOI.VectorOfVariables, SumOfSquares.PositiveSemidefinite2x2ConeTriangle}, Vector{Union{MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.Nonnegatives}, MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle}, MOI.ConstraintIndex{MOI.VectorOfVariables, SumOfSquares.EmptyCone}, MOI.ConstraintIndex{MOI.VectorOfVariables, SumOfSquares.PositiveSemidefinite2x2ConeTriangle}}}}, Union{SumOfSquares.EmptyCone, SumOfSquares.PositiveSemidefinite2x2ConeTriangle, MOI.Nonnegatives, MOI.PositiveSemidefiniteConeTriangle}, MOI.PositiveSemidefiniteConeTriangle, MonomialBasis, MonomialBasis, SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, SumOfSquares.Certificate.NewtonFilter{SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}}}}, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}
| may introduce:
| * Unsupported constraint: MOI.VectorAffineFunction{Float64}-in-PolyJuMP.ZeroPolynomialSet{FullSpace, MonomialBasis, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}
| | bridged by:
| | PolyJuMP.Bridges.Constraint.ZeroPolynomialBridge{Float64, MOI.VectorAffineFunction{Float64}, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}
| | may introduce:
| | * Supported constraint: MOI.VectorAffineFunction{Float64}-in-MOI.Zeros
| * Unsupported variable: SumOfSquares.EmptyCone
| | adding as constraint:
| | * Supported variable: MOI.Reals
| | * Unsupported constraint: MOI.VectorOfVariables-in-SumOfSquares.EmptyCone
| | | bridged by:
| | | SumOfSquares.Bridges.Constraint.EmptyBridge{Float64, MOI.VectorOfVariables}
| | | may introduce:
| * Supported variable: MOI.Nonnegatives
| * Unsupported variable: SumOfSquares.PositiveSemidefinite2x2ConeTriangle
| | bridged by:
| | SumOfSquares.Bridges.Variable.PositiveSemidefinite2x2Bridge{Float64}
| | may introduce:
| | * Supported variable: MOI.RotatedSecondOrderCone
| * Supported variable: MOI.PositiveSemidefiniteConeTriangleIn more details
Consider a polynomial
\[p(x) = \sum_{\alpha} p_\alpha x^\alpha,\]
a vector of monomials b(x) and the set
\[\mathcal{A}_\alpha = \{\,(\beta, \gamma) \in b(x)^2 \mid x^\beta x^\gamma = x^\alpha\,\}\]
The constraint encoding the existence of a PSD matrix Q such that p(x) = b(x)' * Q * b(x) can be written in standard conic form as follows:
\[\begin{aligned} \langle \sum_{(\beta, \gamma) \in \mathcal{A}_\alpha} e_\beta e_\gamma^\top, Q \rangle & = p_\alpha, \quad\forall \alpha\\ Q & \succeq 0 \end{aligned}\]
Given an arbitrary choice of elements in each set $\mathcal{A}_\alpha$: $(\beta_\alpha, \gamma_\alpha) \in \mathcal{A}_\alpha$. It can also equivalently be written in the geometric conic form as follows:
\[\begin{aligned} p_\alpha e_{\beta_\alpha} e_{\gamma_\alpha}^\top + \sum_{(\beta, \gamma) \in \mathcal{A}_\alpha \setminus (\beta_\alpha, \gamma_\alpha)} y_{\beta,\gamma} (e_\beta e_\gamma - e_{\beta_\alpha} e_{\gamma_\alpha}^\top)^\top & \succeq 0\\ y_{\beta,\gamma} & \text{ free} \end{aligned}\]
Should I dualize or not ?
Let's study the evolution of the dimensions m and n of the semidefinite program in two extreme examples and then try to extrapolate from these.
Univariate case
Suppose p is a univariate polynomial of degree $2d$. Then n will be equal to d(d + 1)/2 for both the standard and geometric conic forms. On the other hand, m will be equal to 2d + 1 for the standard conic form and d(d + 1) / 2 - (2d + 1) for the geometric form case. So m grows linearly for the kernel form but quadratically for the image form!
Quadratic case
Suppose p is a quadratic form of d variables. Then n will be equal to d for both the standard and geometric conic forms. On the other hand, m will be equal to d(d + 1)/2 for the standard conic form and 0 for the geometric form case. So m grows quadratically for the kernel form but is zero for the image form!
In general
In general, if $s_d$ is the dimension of the space of polynomials of degree d then $m = s_{2d}$ for the kernel form and $m = s_{d}(s_{d} + 1)/2 - s_{2d}$ for the image form. As a rule of thumb, the kernel form will have a smaller m if p has a low number of variables and low degree and vice versa. Of course, you can always try with and without Dualization and see which one works best.
This page was generated using Literate.jl.